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Abstract 


Recently, there has been increased interest in the possibility of asteroids and comets impacting 
the Earth. These studies have ranged from orbit determination to mitigation, mostly 
concentrating on near Earth asteroids, or NEAs. NEA orbits can usually be determined very 
accurately because many observations are available over several orbital periods. It is possible 
that an asteroid or comet will be discovered on the inbound leg of an impacting trajectory, with 
only a limited amount of time available to collect observations to determine the orbit, and is a 
likely scenario for a long-period comet. Also, unlike NEAs, comet trajectories are perturbed by 
outgassing forces, which are exerted when the comet surface materials are sublimated by the 
radiation of the Sun. The results presented in this paper show the level of orbit determination 
accuracy obtainable for long-period comets discovered approximately one year before collision, 
including non-gravitational perturbations due to outgassing. Comet orbits are designed to impact 
the center of the Earth, and preliminary orbits are determined from simulated observations using 
Gauss’ method. Additional measurements are incorporated to improve the orbit solution through 
the use of a Kalman filter. A comparison is made between observatories in several different 
circular heliocentric orbits, and shows that observatories in orbits with radii less than 1 AU result 
in increased orbit determination accuracy for short tracking durations due to increased parallax 
per unit time. However, an observatory in a circular heliocentric orbit at 1 AU will perform just 
as well if the tracking duration is increased, and the orbit determination accuracy is significantly 
improved if additional observatories are positioned at the Sun-Earth Lagrange points L 3 , L 4 , or 
L 5 . A single observatory at 1 AU capable of both optical and range measurements yields the 
highest orbit determination accuracy in the shortest amount of time when compared to other 
systems of observatories in circular heliocentric orbits at 1 AU. 
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1 Introduction 


The “impact” theory is one explanation of the mass extinctions that have occurred throughout the 
history of the Earth. Many scientists agree that the massive impact crater discovered in the 
Yucatan peninsula is the result of an extraterrestrial collision responsible for the extinction of the 
dinosaurs some 65 million years ago. On a clear night the Moon reveals evidence of the 
devastation caused by asteroid and comet bombardment spanning billions of years, which no 
doubt still occurs today even if it is on a much more infrequent scale. Much media attention has 
been given in recent years to the possibility of an Earth impact, and the attention is not 
unwarranted. The surface of the Earth is still speckled with the scars of past encounters. One of 
the most visible, although certainly not the largest, is the Barringer Meteor Crater in the Arizona 
desert. An object 25 - 60 meters in diameter carved out this one kilometer wide crater upon 
impact 50,000 years ago. In 1908, an object believed to be an asteroid exploded over Tunguska, 
Siberia, destroying an area of forest larger than the city of Los Angeles. It is estimated that an 
impact of this magnitude occurs every one hundred years, and that a globally catastrophic impact 
happens every one million years [1], Existing craters are evidence that our planet has been hit in 
the past. It most definitely will be hit in the future. Nothing could have proved this point more 
than the impact event of comet Shoemaker-Levy 9 with Jupiter in 1994, which left several 
blemishes on the surface of the planet 2 to 3 times the size of the Earth. This event brought to 
our attention a threat that perhaps had not been seriously considered - comets. 

Most current research is in the area of improving orbit determination for known near- 
Earth asteroids (NEAs) with Earth-crossing or Earth-approaching orbits, to predict possible 
collisions many decades into the future. NASA has recently adopted the goal of finding 90% of 
the asteroids 1 km and larger by the year 2008 [2]. This goal may very well be realized, since it 
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is possible to observe the asteroids over many orbital periods due to their close proximity to 
Earth. With the possibility of hundreds of observations over more than one orbital period, the 
orbits can be calculated with sufficient accuracy to predict collision potential decades into the 
future [2]. However, there is no current effort to discover long-period comets that may be 
traveling on Earth-impacting trajectories. A long-period comet is defined as a comet with an 
orbital period greater than 200 years, and could be thousands or even millions of years. The 
orbit determination problem for the long-period comet case is much different than that of 
asteroids, since it is likely that an impacting comet will only be observable for less than half of 
its orbital period before the collision occurs. Astronomers discovered comet Shoemaker- Levy 9 
only 16 months before it collided with Jupiter. It was just enough time to turn the world’s 
telescopes in the direction of our largest neighbor and watch, as the remnants of the comet that 
the massive gravity of Jupiter tore apart hit the planet in a barrage of explosions. The impacts 
left behind distortions in the Jovian atmosphere that were visible from Earth for up to a month. 
Orbit determination for comet Shoemaker-Levy 9 resulted in predictions within 1 1 minutes or 
less of actual impact times for all fragments [3]. 

Orbit determination is the process of collecting a set of measurements and calculating the 
position and velocity of an object at a particular time, or determining a set of six orbital elements 
that define the size, shape, and orientation of the orbit of the object. Preliminary orbit 
determination is performed by using the minimum number of observations required to calculate 
the orbital parameters, and assuming two-body orbital motion about the Sun. The orbit is then 
refined by using methods that can incorporate additional measurements and includes a more 
complete force model including non-gravitational perturbations. The purpose of the current 
study is to investigate the accuracy of orbit determination for comets discovered approximately 
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one year before a collision with Earth occurs. Orbits are constructed so that a collision will take 
place on the inbound leg of the orbit. Preliminary orbits are determined from the minimum 
number of optical observations of celestial latitude and longitude required to calculate position 
and velocity at some epoch. More precise orbits are then determined by adding additional 
optical measurements and observatories, by including radar or laser measurements of range, and 
by incorporating non-gravitational forces due to comet outgassing into the dynamics model. The 
effect of the geometry between the observatory and the comet is also investigated by placing the 
observatory in various locations in circular heliocentric orbits with radii equal to 0.39, 0.72, 1.0, 
and 1.5 astronomical units (AU), corresponding to Mercury, Venus, Earth, and Mars orbits 
respectively. 

This investigation requires a thorough understanding of gravitational attraction and 
orbital mechanics, so this topic is discussed in the next section. Orbit geometries are explained 
and orbital elements are defined in the Appendix. Sec. 3 focuses on warning time and collision 
design criteria. Orbit determination and accuracy are discussed in Sec. 4. Orbit determination 
results obtained by optical measurements of angular position are examined for a family of 
representative comet orbits in Sec. 5, along with results obtained by including additional radar 
information and perturbative effects. Conclusions and suggestions for further study are 
presented in Secs. 6 and 7. 
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2 Orbital Motion 


The iV-body problem in orbital mechanics is the problem of solving the equations of motion for 
N bodies in space, and is formulated using Newton’s laws. The two-body problem is a 
simplification of the A-body problem, for which a complete analytic solution can be found. 
There are numerous approaches to the development of and the solution for the two-body 
equations of motion. The method discussed in the next two sections follows the general 
approach described in detail in Refs. [4] and [5]. It will be seen that the two-body simplification 
requires external forces such as comet outgassing to be neglected, so a method for including 
these perturbative forces in comet orbit determination is described in Sec. 2.4. 


2.1 The A/-Body Problem 

Newton’s laws of motion are [5]: 

First Law Every body continues in its state of rest, or of uniform motion in a straight 
line, unless it is compelled to change that state by a force impressed upon 
it. 

Second Law The rate of change of momentum is proportional to the force impressed, 
and is in the same direction as that force. 

Third Law To every action there is always opposed an equal reaction. 

These three laws of motion were first introduced in 1687 with the publication of Principici, 
although Newton had worked out the equations some 20 years earlier. Also introduced was his 
law of universal gravitation, which states that any two bodies attract one another with a force 
proportional to the product of their masses and inversely proportional to the square of the 
distance between them. Newton formed these laws with the assumption that the bodies were 
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particles, or point masses. Also, the laws are only valid in an inertial reference frame, or a 
reference frame that is at rest or moving with constant velocity. In reality, an inertial reference 
frame does not exist. For practical applications, a reference frame is considered to be inertial as 
long as the observed motion of the particles under study agrees with Newton’s theoretical 
prediction of the motion described by his laws. 

The motion of a system of A bodies can be described analytically using Newton’s second 
law of motion and the law of universal gravitation under the assumptions stated above. The 
gravitational force acting on a body with mass m, due to other bodies in the system with masses 
mi, m 2 , •••, niN is expressed mathematically using Newton’s law of universal gravitation as 


"mm., \ 

V(r-r). i = U.,N, 

j = 1 r n 

j*' 


( 2 . 1 ) 


II 2 2 

where G is the universal gravitational constant and has a value of 6.672x10" Nm"/kg“, r, is the 
position vector from the origin of the inertial coordinate frame to the body on which the force is 
being exerted, r, is the position vector from the origin of the inertial coordinate frame to the body 
exerting the gravitational force, and r,y is the distance between the i-th and / - th bodies, where 

( 2 - 2 ) 


Newton’s second law of motion states that the total force acting on the body is 


f = m 


d 2 r ; 


! ' 


dr 


m , 


dVj 

dt 


(2.3) 


d r 

where — 7 - and v, are the inertial acceleration and velocity of the body, respectively. 

dr 
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Assuming no other external forces are acting on the system, the expression describing the motion 


of a system of N particles is given by 



(2.4) 


There is no known general solution to Eq. (2.4). A complete solution requires 6 N integrals, of 
which only 10 are obtainable. The 10 constants of integration result from the fact that total linear 
momentum, angular momentum, and energy of the system are conserved. From the conservation 
of total linear momentum it can be shown that the center of mass (CM) of the system moves with 
constant linear velocity, and that the position of the CM as a function of time is 

r{t)=\ Q {t-t Q )+r 0 , (2.5) 

where Vo and ro represent 6 of the constants of integration. The total angular momentum, h, of 
the system is 

N 

h = X r / Xm /V / , (2.6) 

i = 1 

which can be shown to be conserved so that h represents three more of the constants of 

integration. The tenth constant is the total system energy, t, which is the sum of the total kinetic 

and potential energy of the system. Complete derivations and proofs of the conservation of total 
linear momentum, angular momentum, and system energy are presented in detail in Ref. [4] . 

An analytical solution can be found for Eq. (2.4) for the relative motion of two particles 
under their mutual gravitational attraction. All other forces acting on the particles, including 
iV-body gravitational attraction, can be described in terms of perturbations. For this reason, the 
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simplified approach of solving the equations of motion for two bodies is ideal for describing the 


motion of comets in orbit about the Sun. 


2.2 The Two-body Problem 

The inertial reference frame containing the two bodies (the Sun and the comet) is defined by 
three orthogonal unit vectors, i x , i , and /, . Positive i is normal to the ecliptic plane and 

points in a direction parallel to the orbital angular momentum vector of the Earth, positive i x 

points in the direction of Earth vernal equinox and lies in the ecliptic plane, and positive i lies 

in the ecliptic plane and completes the orthogonal coordinate system. The mass of the comet is 
negligible compared to the mass of the Sun. The only perturbative effects to be included in 
determining the comet orbits are those due to outgassing because they are unknown; 
gravitational attraction from the planets in the solar system, including Earth, can be modeled 
very precisely and therefore are not considered. 


/x 




Figure 1. Inertial reference frame unit vectors 
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The equations of motion for two bodies are given by Eq. (2.4) with N = 2 and i = 1,2, 


i.e., 


and 


m . 


d 2 r l 

dt 2 


r n 


777 , 


d r 2 

Jt 2 


Gi 


-r t ) 

-r 2 ). 


(2.7) 


(2.8) 


The equation of relative motion for the two-body problem is found by subtracting Eq. (2.8) from 
Eq. (2.7) to obtain 



= -G 



2 ) 


r , 


(2.9) 


where r = r 2 - r, , and r is the magnitude of r. The term G{m x + m 2 ) is called the gravitational 
parameter and is represented by //. With this substitution, Eq. (2.9) becomes 



+ 4r = 0. 
r~ 


( 2 . 10 ) 


Eq. (2.10) is a second order, nonlinear, vector differential equation describing the relative motion 
of two bodies. Even though it is nonlinear, it is possible to find a completely analytic solution to 
this equation through six independent constants of integration. These six constants can be the 
components of position and velocity, or they can be a set of six parameters that define the orbit 
of the object around the primary body, called orbital elements. Orbital elements are described in 
detail in the Appendix. 
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2.3 Orbit Propagation 


If the position and velocity of an object are known at some time to, the position and velocity can 
be calculated at any time t using the expressions 


r = /r 0 + gv 0 , 

(2.11) 

v = /r 0 + gv 0 , 

(2.12) 


where the variables / and g are called Lagrange coefficients, also known as the /and g functions. 
There are many ways to calculate the Lagrange coefficients, one of them being through the use 
of universal functions, which are used extensively throughout Ref. [4] and summarized in the 
equations to follow. The Lagrange coefficients in terms of universal functions are 


where 


/= 


8=~T={ r 0 U l +(7 o U 2 )> 

v/A 

rr n 


■ , U 2 

8 = 1 

r 


<?o =^( r o- v o)- 


(2.13) 

(2.14) 

(2.15) 

(2.16) 

(2.17) 


The universal functions are expressed differently whether the orbit is parabolic, hyperbolic, or 
elliptic. The expressions for the first four universal functions for elliptical orbits are 
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U 0 = cos (4zx) , 


(2.18) 


U i 


sin (4z%) 


U 2 = 


1-cos (4zZ) 


U, = 


yfzZ - sin (Vz/jf) 

z4z 


(2.19) 

(2.20) 
(2.21) 


where z is the reciprocal of the semi-major axis of the orbit and can be expressed as 


z 


M 


(2.22) 


The variable X ' s an independent variable used instead of time, sometimes referred to as the 
generalized anomaly, and is related to the eccentric anomaly by 


X = 



(2.23) 


The universal functions are used to obtain the generalized form of Kepler’s equation 


y[M( t ~ t o)~ r o^i +(7 o^2 + U 3 - (2.24) 

If time is given, Eq. (2.24) must be solved for x to find the position in the orbit. This equation is 
transcendental in X-> therefore an iterative scheme is necessary. Ref. [5] suggests using a first 
guess of 

X g =4 m z(t~t 0 ), (2.25) 
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where t is the time at which r and v are to be calculated and to = 0. This guess is used in a 


Newton iteration algorithm 


t-t 


Z g+ i =Z g + 


dt 


dX 

x=z g 


(2.26) 


where % g is first used to calculate the universal functions and then Eq. (2.24) is used to calculate 

t g by setting t = t g . Substituting universal functions into Eq. (4.4-17) of Ref. [5] yields an 

. , dt . 

expression for the term — , i.e., 

dX 


dt 

dX 



(r 0 U 0 + (T 0 U X + U 2 ) . 


Eq. (2.26) is iterated until t - t g converges to the pre-assigned value. 


(2.27) 


2.4 Perturbing Forces 

The solution to the two-body problem presented in the previous sections describes the relative 
motion between two particles when the only force acting on the particles is their mutual 
gravitation. In reality, other forces act on the bodies and perturb the motion. Some examples of 
these perturbing forces are gravitational forces from other bodies in the system or the asphericity 
of the central body, atmospheric drag, electromagnetic forces due to interaction with magnetic 
fields, or solar radiation pressure. The additional forces and the deviations from two-body 
motion are called perturbations, and there are two methods for analyzing them. In the method of 
special perturbations, the equations of motion, including all perturbing accelerations, are 
numerically integrated. In the method of general perturbations, an analytical solution is achieved 
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through variation of parameters. Both methods are described in Ref. [5]. The method of special 
perturbations will be used to include the perturbative effects due to comet outgassing on the two- 
body comet trajectories under study. 

A comet is often thought of as a fuzzy point of light in the night sky. A comet appears 
“fuzzy” because the heat radiated from the Sun is sublimating the comet surface, and the 
evaporated materials form a cloud around the core, or nucleus, of the comet. The cloud of 
evaporated materials is called the coma. Solar radiation pressure and the solar wind blow this 
dust away from the nucleus and a tail is formed. The composition of the nucleus is still mostly 
unknown, but it is generally thought of as a “dirty snowball”. Based on spectroscopic 
observations of comet tails it has been speculated that the nucleus is made mostly of frozen water 
and organic and silicate compounds [6]. When the comet gets close enough to the Sun so that 
the materials in the nucleus begin to evaporate, the comet is said to be outgassing. This 
outgassing exerts forces on the comet that perturb the orbit, so it is necessary to determine 
whether the magnitude of these deviations will make it impossible to determine the orbit to some 
specified level of accuracy. 

A commonly used method for modeling outgassing accelerations was developed by 
Marsden, Sekanina, and Yeomans, and is described in Ref. [7]. By carefully studying the orbits 
of over 20 periodic comets they were able to derive the following expression for acceleration due 
to outgassing, 

a = g(r)[/f ,r + A 2 t + /4 3 n], (2.28) 

where 




f \ 

—m 


( \ 

n 

r 


1 + 

f r 


— 


— 


UJ 



l r o J 



—k 


(2.29) 
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The unit vector r points in a direction parallel to the Sun-comet vector r in the orbit plane, t is 
perpendicular to r in the orbit plane, and n is perpendicular to the orbit plane. This model is 
based on the assumption that the comet is an “icy conglomerate” - an object consisting mostly of 
frozen water that holds together bits and pieces of rock-like material. This model also was used 
to generate accurate ephemerides for Halley’s comet [8]. In Eq. (2.29) the variable r is the 
heliocentric distance of the comet in AU, and all other parameters are constants based on studies 
of vaporization rates of comet nuclei materials and are derived in Ref. [7]. The constants are 
reported to be m = 2.15, n = 5.093, and k - 4.6142. The normalizing constant c is defined such 
that g(l) = 1, and r 0 is the scale heliocentric distance of high outgassing activity with a value 
equal to 2.808 AU for frozen water. Substituting these values into Eq. (2.29) and solving for the 
normalizing constant yields c = 0.1113. In Eq. (2.28), A, (i = 1, 2, 3) are constants describing the 
components of acceleration in the radial, transverse, and normal directions. Values of A t and A 2 
calculated by studying changes in comet orbital periods are listed in Table I of Ref. [7] for 
various comets. Values for A3 are not given because the acceleration force normal to the orbit 
plane had no detectable effect on the orbital periods studied. However, a normal acceleration 
force may be present if the comet is rotating or tumbling, and this component must be included 
in simulations designed to predict collisions. Without any given data for such a force it is 
assumed that A3 = A2. It should also be noted that the model is based on comets that have made 
multiple solar passes. A comet approaching the Sun for the first time may be subject to larger, or 
smaller, outgassing accelerations than described by this model, depending on the composition of 
the nucleus. 

Before undertaking orbit determination simulations it is important to understand the 
magnitude of displacement from two-body motion along the comet trajectory caused by 
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outgassing accelerations. These accelerations are included as an additional term in Eq. (2.10) so 


that the equation of motion for an outgassing comet is 


d r _ // 

dt 2 ~ r 3 


r + g(r)kr + A 2 f + A 3 h 


(2.30) 


By comparing results obtained through numerical integration of Eqs. (2.10) and (2.30), the 
perturbation to two-body motion caused by comet outgassing is determined as 


fr{t) = r{t) n -r{t) o . 


(2.31) 


where r„ is the position of the comet calculated through integration of Eq. (2.10), and r„ is the 
position calculated by integration of Eq. (2.30). The integration is performed for a sample long- 
period comet (orbital period greater than 200 years) on an Earth-impacting trajectory with orbital 
elements r a = 300 AU, r p = 0.1 AU, i = 45°, Cl = 30°, CO- 143.2°, and v Q = -163.9° (criteria for 
creating such orbital elements is presented in Sec. 3). The orbital period for this comet is 1839.4 
years. The equations of motion are integrated from to to tf, corresponding to times when the 
comet is at heliocentric distances of 5 AU and 1 AU, respectively. The maximum values for the 
acceleration coefficients listed in Ref. [7] are used for this analysis. These values are Ai = 3.61 
AU/10 8 day 2 (5.4 km/day 2 ) and A 2 = A 3 = 0.3269 AU/10 8 day 2 (0.49 km/day 2 ). Fig. 2 shows the 
absolute value of each component of Sir, denoted by Sr x , Sr y , and Sr z , versus time on a 
logarithmic scale. Also shown is the magnitude of Sr versus time, Sr. which is approximately 
2,500 km by the time the comet reaches a heliocentric distance of 1 AU for this sample case. 
This deviation is sufficiently large compared to the size of the Earth that outgassing must be 
included in orbit determination accuracy calculations. However, the more important result is that 
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the perturbations are less than 10 km until about 200 days, and can be considered negligible. At 
200 days the comet has traveled to a heliocentric distance of about 2.5 AU. Since these 
deviations are so small compared to the heliocentric distance of the comet, a first order 
perturbation adequately represents the motion of the comet, and a linear model is adequate to 
represent the perturbations for orbit determination purposes. This is confirmed through Monte 
Carlo simulations, the results of which are presented in Sec. 4.2.4. A similar result is obtained 
for a short-period comet (orbital period less than 200 years) with orbital elements r a = 40 AU, 
r p = 0.1 AU, i = 50°, £2 = 30°, (O = 143.6°, and Vo = -164.8°. The orbital period of this comet is 
89.8 years. Fig. 3 also shows negligible perturbations until approximately 200 days, when the 
comet has reached a heliocentric distance of 2.6 AU. 






Figure 2. Deviation in position due to outgassing forces for a long-period comet 
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Figure 3. Deviation in position due to outgassing forces for a short-period comet 


3 Impact Trajectory Criteria 


The accuracy of orbit determination depends on many factors. Some of these that are important 
to this investigation are the number of observations available, the time interval over which those 
observations are made, the resolution of the measuring device, and the geometry between the 
observatory and the object. A standard method for examining these effects is to construct a 
collision orbit, use this “reference” orbit to calculate a set of perfect measurements and add 
random noise to them, and then determine the “estimated” orbit using the noisy measurements 
with the orbit determination method described in Sec. 4. An important consideration in 
designing these true orbits is warning time. Warning time is the interval between the time when 
the comet orbit becomes known and the time of collision, and is greatly dependent on the 
heliocentric distance of the comet at the time observations are made. Warning time is discussed 
in Sec. 3.1, with the purpose of establishing an approximate distance at which the comet should 
be discovered to supply enough time to observe and determine the orbit to the desired degree of 
accuracy, while leaving roughly one year to plan and execute an evacuation scenario or even 
some type of mitigation technique. A set of comet orbits representing several different orbit 
geometries are designed to guarantee Earth-impact using the equations presented in Sec. 3.2 [9]. 


3.1 Warning Time 

As stated previously, warning time is the interval between the time when the comet orbit 
becomes known with sufficient accuracy and the time when the collision takes place. It is 
possible that the comet will travel through several orbits in this time, but the possibility also 
exists that the comet will be discovered on its terminal orbit, greatly reducing the amount of 
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warning time. A system designed to detect impacting comets should be able to provide enough 
information to determine such an orbit, and it is important to calculate the amount of warning 
time to be expected in this case. 

Fig. 4 shows the Earth in a circular heliocentric orbit of radius r c - 1 AU, and a comet in 

A A 

a coplanar elliptical orbit which crosses the orbit of the Earth. The i x and i axes lie in the 

ecliptic plane, with positive i x in the direction of vernal equinox. The comet is discovered after 

aphelion, and the collision takes place before perihelion. The time between discovery and 
collision is the quantity to be determined. The motion of both the comet and the Earth about the 
Sun is governed by the two-body orbit mechanics described in Sec. 2.2. 
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Inclination does not affect orbital period or the time of flight between one point and 
another; therefore, studying the coplanar case is suitable for generalizing warning time analysis. 
Perturbative effects on the comet due to the gravity of the Earth and comet outgassing are not 
considered in the analysis to follow, since these effects will not significantly change the time of 
flight on the relatively short time interval being investigated (compared to the orbital period of 
the comet). 

The eccentric anomaly, E , at any point on an elliptic orbit is given by Eq. (A14). If E 
denotes the eccentric anomaly at the point of discovery where r is the distance between the Sun 
and the comet, and E c denotes eccentric anomaly at the point of collision where r = r c = 1 AU, 
then the time of flight between these two points is calculated using Eq. (4.2-9) of Ref. [5] given 
here as 

t-t c = ^j— [(£-csin £)-(£). -esin £) )]. (3.1) 

Eq. (3.1) is used to calculate the warning times for a sample family of short- and long- 
period comets whose orbits become known at distances of 5, 6, and 7 AU. The family is 
constructed with perihelion distances, r p , equal to 0.1 and 1 AU, and aphelion distances, r a , equal 
to 15, 20, 25, ..., 100, 200, 300, ..., 1000, 2000, ..., 50xl0 3 AU. The range of orbital periods for 
these values of aphelia is between 20 years and 4xl0 6 years, reaching out to the middle of the 
Oort cloud. Long-period comets are typically defined as those with periods greater than 200 
years. 

Warning time as a function of r a is displayed in Fig. 5. Warning time does not change 
significantly for aphelia greater than 1000 AU, and reducing r p by a factor of 10 results in a loss 
of warning time by about the same amount as reducing the detection distance by 1 AU. The 
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Figure 5. Times of flight between detection and collision for short- and long-period comets 


range in warning times is between a little more than 2 years for the case where the comet is 
discovered at a distance of 7 AU, with r p equal to 1 AU, and 9.5 months for the case where the 
comet is discovered at a distance of 5 AU, with r p equal to 0.1 AU. 


3.2 Designing Earth-Impacting Orbits 

By design, all of the orbits to be studied pass through the ecliptic plane at a heliocentric distance 
of 1 AU. If the ecliptic plane is the reference plane, a comet in an inclined orbit is at the 
ascending or descending node when passing through the ecliptic. Therefore, the collision takes 
place when the comet is at the ascending or descending node, and the heliocentric distance r is 
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equal to r c . Referring to Fig. Al, we see that v = -co when the object is located at ascending 
node, and v= n-CO when the object is located at descending node. The equation of the orbit can 
be written in the form 


r p { 1 + g) 
1 + ecosv 


(3.2) 


which after substituting the collision requirements becomes 

r (l + e) 

y — — ? 

l±<?COS CO 


(3.3) 


where the positive and negative signs are associated with the ascending and descending nodes, 
respectively. Solving for argument of periapsis, CO, gives 


± COS CO- — 

r. 


Ail 

V e) 


1 

e 


(3.4) 


The comet orbit will intersect the orbit of the Earth as long as 


and 


< 


r < r . 

c — a 


(3.5) 

(3.6) 


Therefore, Eq. (3.4) can be used to solve for the value of ft) that ensures a collision to take place 
at the ascending or descending node, as long as Eqs. (3.5) and (3.6) are satisfied. When the 
comet orbit lies in the ecliptic plane, the ascending and descending nodes will be undefined since 
i = 0°. In this case, the value of co is arbitrary because a collision is ensured to occur as long as 
the conditions in Eqs. (3.5) and (3.6) are satisfied. 
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4 Orbit Determination and Collision Probability 


It was shown in Sec. 2.2 that six parameters are necessary to define an orbit, either a set of six 
orbital elements or the six components of position and velocity. Therefore, a minimum set of six 
measurements is required to calculate an orbit, usually two angular measurements at three 
separate observation times. Newton was the first to describe a method for finding the orbit of a 
body from three observations, and Edmund Halley was the first person to apply this technique to 
the calculation of orbits. He used it to prove that a comet observed by Kepler in 1607 was the 
same comet he observed in 1682, and predicted that it would return in 1758. On Christmas day, 
1758, the comet appeared in the night sky, fulfilling his prediction and ever after bearing his 
name [5]. 

Preliminary orbit determination is the process of using the minimum number of 
observations to calculate a first approximation of the six orbit parameters required to define an 
orbit. This approximation can then be combined with additional angular measurements, range 
information, and effects due to perturbative forces to obtain the orbit more precisely. The comet 
orbits under study have been designed to cross the ecliptic plane at a heliocentric distance of 
1 AU where a collision at the center of the Earth is assumed to occur. The designed orbit defines 
a nominal set of observations consisting of celestial longitude, (p, celestial latitude, A, and in 
some cases range, p , to which random noise is added to simulate observations collected from an 
observatory. The accuracy of the orbit determination can be quantified by comparing the orbits 
calculated with noisy measurements to the nominal orbit for different tracking schedules and 
observatory locations. Orbit determination methods are discussed in Sec. 4.1, and orbit 
determination accuracy is discussed in Sec. 4.2. 
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4.1 Methods of Orbit Determination 


Ultimately it would be Gauss who would be remembered for solving the orbit determination 
problem based on three observations, although he knew himself that his method was only an 
approximation to be used for preliminary orbit determination. He’s quoted as stating that 
improvement could only be made by accumulating “the greatest number of the most perfect 
observations, and to adjust the elements, not so as to satisfy this or that set of observations with 
absolute exactness, but so as to agree with all in the best possible manner” [5]. His method of 
weighted least squares is one way of combining a “great” number of observations to calculate an 
improved estimate of the orbit. This method usually employs observational data collected 
together as a “batch”, and is applied to minimize the sum of the squares of the weighted residuals 
of the measurements. The batch process requires all data to be stored, so when a new 
measurement becomes available, it is added to the batch and the entire set is either reprocessed 
or augmented using a sequential batch process. Simple formulas exist for adding a new 
measurement without having to reprocess the entire data set, and various levels of importance 
can be placed on the measurements through weighting factors. However, the method of weighted 
least squares has no method for incorporating system noise to account for incomplete force 
modeling or measurement noise to account for imperfect observations. An orbit determination 
method that does incorporate system and measurement noise statistics makes use of the Kalman 
filter, a sequential data processing algorithm that yields a minimum error variance, or minimum 
trace of an error covariance matrix. This method is described in Sec. 4.1.2. 
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4.1.1 Preliminary Orbit Determination Using Gauss’ Method 

Gauss’ method makes use of measurements of longitude, cp , measured from vernal equinox in 

the ecliptic plane, and latitude, A , measured from the ecliptic plane toward the north celestial 
pole. The tilde superscript denotes the fact that the measurements are assumed to have some 
error associated with them due to the measurement resolution of the observing device. 
Following the algorithm presented in Sec. 5.8 of Ref. [5], assuming that optical measurements (p 

and A are available at three times t\, h, and p, a unit vector along the line of sight between the 
observatory and the comet is expressed as 

L,. = (cos (p j cos A t ) i x + (sin (p j cos A, ) i + (sin \ ) i , , i = 1, 2, 3 . (4.1) 

The heliocentric position of the comet can be expressed as 

r = d + L/?, (4.2) 

where d is the heliocentric position of the observatory, and p is the magnitude of the distance 
between the observatory and the comet. If the vector r(t) is expressed in terms of the Lagrange 
coefficients /and g at h. then 

r(0=/(i2.v 2 ,f-f 2 ) r 2 +g(r 2 ,v 2 ,f-f 2 ) v 2 , (4.3) 

and with 

/; =/(r 2J v 2 ,r f ~t 2 ) (4.4) 

Si =g(r 2 > y 2 > t i -h)’ ( 4 - 5 ) 
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a set of nine equations in the nine unknowns r 2 , v 2 , p\, pi, and p can be expressed as 


/i r 2 + <?1 V 2 = d l +A L 1. 

(4.6) 

**2 — d 2 T PA 2 > 

(4.7) 

fPi+SPi =d 3 +A L 3- 

(4.8) 


This set of equations is reduced to six independent equations in the six unknowns r 2 and v 2 by 
cross multiplying the i-th equation with L, and eliminating the i, components to obtain the 
following set of equations written in matrix notation, 
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(4.9) 


where r x , r y , and r z are the components of r 2 , and v x , v y , and v z are the components of v 2 . The 
variables / and g are functions of r 2 and v 2 and are obtained through the process outlined in 
Sec. 2.3. Eqs. (4.9) are nonlinear, and the procedure for solving this set of equations is: 

1. Estimate the vector r 2 using Eq. (4.2) at f 2 , where d 2 and L 2 are known, and p 2 is 
assigned a value of 5 AU (the distance at which the comet is assumed to be discovered); 

2. Estimate the vector v 2 by forming the cross product of Li and L 3 to obtain a unit vector 
along the direction of the angular momentum, h ; form the cross product of h with a unit 
vector in the direction of the estimated position, r 2 , to obtain a unit vector in the direction 
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of the estimated velocity, v 2 ; estimate the magnitude of v 2 using the expression for 
circular velocity, i.e., 



3. Using these estimates, compute the values of f\, gu / 2 , and g 2 using the equations 
presented in Sec. 2.3, and substitute into Eq. (4.9) to solve for the components of r 2 and 

v 2 ; 

4. Update /i, gufi, and g 2 and iterate until r 2 and v 2 converge. 

In all of the analysis to follow, Gauss’ method is used with three observations spaced seven 
days apart to determine preliminary orbits and other required initial parameters for the Kalman 
filter. The preliminary orbit serves as an initial guess for the filter, and additional observations 
are collected to attain more accurate orbit determination. 

4.1.2 Final Orbit Determination Using the Kalman Filter 

As previously stated, the Kalman filter is a sequential filter, or a recursive data processing 
algorithm in which an estimate of the state, namely the six components of position and velocity, 
determined from measurement data available at the current time is updated when measurements 
become available at subsequent times. The objective of any type of filter is to obtain this 
estimate while minimizing the error in some respect. The Kalman filter is an optimal estimator 
that minimizes the error variance, or the trace of the error covariance matrix. The Kalman filter 
can incorporate all available measurements and measurement uncertainties, along with the 
dynamics models of the comet and observatory and the uncertainty in these models to produce an 
estimate of comet position and velocity at the time of each observation. 
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4.1. 2.1 State Propagation 


Using a Kalman filter requires some basic assumptions to be made about the problem. First, the 
problem must be described using a linear model. Since the orbit determination problem is 

5|s 

nonlinear, the approach is to linearize about a reference trajectory. Let r (t) and v (!) be the 
position and velocity of the comet on the reference trajectory. The comet will not follow this 
reference trajectory exactly, and the actual position, r(t), and velocity, \(t), will deviate from 
r ( t ) and v it). Expanding r(t) and v(/j in a Taylor series about the reference trajectory yields 
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where [ ] indicates that the matrix of partial derivatives is evaluated on the reference trajectory, 
higher order terms are neglected, and 


*o = r{t 0 )-r'{t 0 ), 

<fr 0 = v (h>)- v *(0- 


If a state deviation vector x is defined as 


it)= 


Mt) 

Mt) 


then Eqs. (4.10) and (4.11) can be written as 


x{t)=<f>{t,t 0 )x{t 0 )+w{t 0 ). 


(4.12) 

(4.13) 


(4.14) 


(4.15) 
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where <t>(/,/ 0 ) is the state transition matrix, given by 




M<) 

MM 
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(4.16) 


As before, [ ] indicates that the matrix of partial derivatives is evaluated on the reference 
trajectory. The numerical values of the elements populating <f> can be calculated from the 
analytical expressions of Eqs. (9.84) through (9.87) in Ref. [4], which are repeated below, 
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= — (v - v 0 )(v - v () f +~'j k (l - / ) rr ( ; -CrVo]+gI, 

(4.20) 


where C is introduced for convenience of notation and is defined in terms of universal variables 
by 


C = 



(3 u 5 -zU 4 )-{t-t 0 )U 2 . 


(4.21) 


The quantity C is determined using the equations given in Sec. 2.3, along with the following 
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expressions for the universal variables U 4 . and U 5 , 


u*=u l u 1 ~(ul-zul). 


(4.22) 


U 5 = 


1 


u,-x 




(4.23) 


When %= U\ =0. it can be shown through identities given in Ref. [4] that U5 = 0. 


4.1. 2.2 The Observation Equation 

The observations are assumed to occur at discrete points in time and follow the linear relation 


y* 


_ 5 


(4.24) 


where the subscript k denotes the discrete point in time, f* . The vector y * is an mx 1 observation 
residual vector, the elements of which are the difference between the measurements on the actual 
trajectory and the measurements on the reference trajectory, and e* are the errors in the 
measurements. The matrix H is an mx 6 matrix of partial derivatives of the observations with 
respect to the state vector, where the state vector is a 6-dimensional column vector of the 
components of position and velocity at time t±. This matrix is generally based on the geometry 
of the problem and is sometimes referred to as the observation matrix. For measurements of 
longitude and latitude, the observation matrix is expressed as 
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where p is the position vector of the comet with respect to the observatory, i.e., 


P = PJx + PyK + pjz 


(4.26) 


and p is the magnitude of p. When each observation consists of two angular measurements, the 
dimension of y* is 2x1 and the dimension of H is 2x6; if a measurement of range is added to 

each observation, the dimension of y* is 3x1 and the dimension of H is 3x6. The first two rows 

of H will be the same as above, and the last row will contain the partial derivatives of the 
measurement of range with respect to the state, i.e., 
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(4.27) 


The second Kalman filter assumption is that the measurement noises are white, meaning 
the values of noise are not correlated in time. In other words, the value of noise at a particular 
time cannot be predicted based on knowledge of the noise value at previous times. Therefore, 
the measurement errors £/ ( in Eq. (4.24) are uncorrelated errors with 

e{ e(f) }= 0 

e{ e(<) e(i) T }=R(»), (4.28) 


where the mxm matrix R(7) is the covariance of £ and describes the uncertainty in the 
observations. 
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4.1. 2.3 Time and Measurement Updates 

The initial state deviation x(to) will not be known precisely a priori. Therefore, the dynamics 
model of Eq. (4.15) is written in terms of an estimate of the initial state deviation x(t 0 ), i.e., 

\{t) = 4>(f , t 0 ) x(t 0 ) + w(f 0 ) , (4.29) 


where w (t 0 ) is white noise with 


£■{ w(t) }= 0 

e { w(r) w(t) T }=Q(t). (4.30) 

The 6x6 matrix Q(t) is the covariance of w and describes the uncertainty in the dynamics model. 
A 6x6 covariance matrix P is defined in terms of the error in the estimate of the state deviation 
vector, 

e { [x(f)-i(r)] [x(r)-x(t)] T }= P . (4.31) 

The Kalman filter is used to estimate and update the state deviation vector from the time 
when one measurement is processed, E-i, to the time the next measurement is processed, Z>. This 
process can be thought of as consisting of two parts: a time update and a measurement update. 
The procedure uses the solution at time fk-i as the starting point for the update at time Z>. 
Consider the state deviation at a time just after a measurement update, denoted by t k _ t , and a 

time just before the next measurement update, t~ k . The time updates of the state deviation vector 
and the covariance matrix P are given by Eqs. (4.2-17) and (4.2-18) in Ref. [10] repeated here, 
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(4.32) 


K = ®{t k ,t k _ l ) r k _ { , 

p; = <t>{t k ,t k _ 1 ) p;_ t <t>{t k ,t k _J +Q h . (4.33) 

The measurement updates of the state deviation vector and the covariance matrix P are given by 
Eqs. (4.2-5) and (4.2-12) in Ref. [10] as 

(4.34) 

P. t =(l-K 1 H 1 )P i :(l-K 1 H,) T +K 1 R,Kl, (4.35) 

where K is the Kalman gain, given by Eq. (4.2-15) in Ref. [10] as, 

K,=P;Hl[H,P,-Hj+R [ ] _, f (4.36) 

and is an optimal gain matrix that weights the measurements more heavily if K is large and less 
heavily if K is small. 

The filter generates the best prediction of the state deviation x( at time f* before the 
measurement at time t k is processed, and generates the best prediction H A x“ of what the 
measurement deviation will be before the measurement is actually taken. The measurement 
residual is calculated as the difference between the measurement deviation y * and H A x A , and 
then multiplied by the gain matrix, K*, to produce a correction term to be added to x A to obtain 

x k . This correction term is proportional to both K* and the measurement residual [y A — H A x A ] . 
The algorithm has a predictor-corrector structure. The Kalman filter algorithm is completely 
defined by specifying <f> and H for all times of interest, and supplying x(f 0 ), y*, and the 
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uncertainties Po, Q it), and R(/) [11], 

An initial estimate of the state deviation x(/ 0 ) is supplied in the form of a 6x1 column 
vector of deviations in position and velocity from the true comet orbit, and a corresponding 
initial 6x6 covariance matrix Po describing the uncertainties in \(t Q ) . The initial conditions are 
obtained from a preliminary orbit determined by Gauss’ method, using three observations 
collected by an observatory with measurement resolution I] equal to 0. 1 arcseconds (comparable 
to Hubble Space Telescope). The Gauss algorithm presented in Sec. 4.1.1 is used 100 times with 
100 sets of six observations that have normally distributed noise with zero mean and standard 
deviation equal to //. The samples are truncated at ±2 tj. From this process, 100 different 
preliminary orbit solutions are obtained. The initial estimate of the state deviation, x(t 0 ) , is the 

difference between the Gauss solution and the designed orbit. The initial covariance matrix is 
calculated from the 100 state deviation vectors. 

Since outgassing forces have been shown to be negligible during the time interval when 
observations are being collected, all time updates are performed by way of the two-body state 
transition matrix of Eq. (4.16). Comet outgassing introduces an uncertainty in the two-body 
assumption and is modeled as two effects - a varying uncertainty that changes with time and a 
constant uncertainty that acts over the entire orbit. The time-varying uncertainty is modeled as 
white state noise with zero mean and diagonal covariance Q(t). The elements of Q it) are 
calculated with the aid of Eq. (2.28). Since it is impossible to know ahead of time the values of 
the outgassing coefficients for a newly discovered comet, a distribution of values between the 
maximum and minimum values of A, (; = 1, 2, 3) given in Ref. [7] is used to determine standard 
deviations for the coefficients, <7 A . (i = 1, 2, 3). Using the o A . (i = 1, 2, 3) in place of the A,- (i = 1, 

2, 3) in Eq. (2.28), outgassing acceleration a is determined. Assuming a is constant over the 


33 



time interval between observations, the deviation in position due to outgassing acceleration is 


expressed as 


Sr = ^a'{t k - ) 2 , 


(4.37) 


and the deviation in velocity is expressed as 


<5v — ot (t k t kl ) , 


(4.38) 


where the quantity (t k - t k , ) is the time interval between observations. The uncertainty in the 
system model can be expressed as 


Q,U = <S<7. i = l, 2. 3 

(4.39) 

-' = 4,5,6. 

(4.40) 


The constant element of uncertainty caused by outgassing is described in Sec. 4.2.2. 


4.2 Orbit Determination Accuracy 

Body plane (B-plane) targeting is a method commonly used for determining spacecraft closest 
approach distance to a target planet during the design phase of interplanetary missions, and for 
interplanetary orbit determination and maneuver targeting. This method is easily applied to 
define a metric to quantify orbit determination accuracy for Earth-impacting comets in terms of 
miss distance. Results from B-plane analysis can be used further to determine Earth collision 
probability, as in Ref. [12]. 
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4.2.1 The B-plane 

The B-plane is defined as a plane passing through the center of the target planet, perpendicular to 
the asymptote of the incoming hyperbolic orbit. The vector B is a vector in that plane, on a line 
from the center of the target planet to the point where the asymptote intersects the plane, as 
shown in Fig. 6. The vector B determines where the point of closest approach would be 
assuming the target planet had no mass and did not influence the trajectory. A coordinate system 

with origin at the center of the target planet is defined by three orthogonal unit vectors T , R , 

a, A, yv 

and S . The unit vector S points in the direction parallel to the incoming asymptote, T is normal 

A, A. 

to S and lies in the ecliptic plane, and R completes the orthogonal basis. The construction of 
the basis is such that S is perpendicular to the B-plane, and T and R lie within it. 
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In the present study, the B-plane is constructed with the Earth as the target planet. 
Ignoring the gravity of the Earth, the velocity at infinity can be calculated as 

v ~ = v c (f c )-v £ (f c ), (4.41) 

where v c (t c ) and v E (t c ) are the velocities of the comet and the Earth in the inertial reference 
frame at the designed time of collision, respectively. Thus, 

S = — . (4.42) 

V oo 


yv yv 

The unit vector T lies in the ecliptic plane perpendicular to S , therefore 



The unit vector R is 


(4.43) 


R = SxT. 


(4.44) 


As stated previously, the vector B is a vector in the B-plane, extending from the center of 
the Earth to the point where the hyperbolic asymptote intersects the B-plane. Since S is normal 
to the B-plane, and since both B and S lie in the trajectory plane, B must satisfy the expressions 


Bx(i 

'-s) 

!=h, 

(4.45) 

B ■ (v 

’-s) 

l= o, 

(4.46) 


where h is the angular momentum of the comet relative to the center of the Earth, 
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h = r oo xv to , 


(4.47) 


and Too is the relative position of the comet with respect to the Earth at the designed time of 
collision, 

E, = r c (t c )-r h -{t c ). (4.48) 


Forming the cross product of S with Eq. (4.45) gives 

Sx(bxs)= — (sxh), 


or 


(s-s)b-(s-b)s = — (sxh). 


(4.49) 


(4.50) 


Since S is a unit vector, and since B and S are mutually orthogonal, Eq. (4.50) becomes 


B = — (sxh). 

v„ 


(4.51) 


The magnitude of B is the miss distance. It is common to describe B using the two components 
along the T and R directions, 


and describe them verbally as “B dot T” and “B dot R”. 


B T — B ■ T , 

(4.52) 

B r = B ■ R , 

(4.53) 

“B dot R”. 
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4.2.2 State Uncertainties at Nominal Collision Time 


A collision occurs at the center of the Earth by design. Therefore, if the collected observations 
were perfect, producing perfect orbit determination, B would be zero. Of course, the 
observations have random noise associated with them and a perfectly determined orbit is not 
expected. The objective of the Kalman filter process is to determine a vector B and an associated 
uncertainty. This uncertainty is determined by propagating the final Kalman covariance matrix 
P t + from t k to t c by 

P t =<%,(,) (4.54) 

with a modification required for the calculation of Q/. The state noise covariance given in 
Eqs. (4.39) and (4.40) is valid for short time intervals over which outgassing forces are assumed 
constant. The interval between the time of the last observation and the designed time of collision 
is too long for this assumption to remain valid (approximately one year), so this long interval is 
divided into 20 shorter intervals (each roughly two weeks). The modified time update is 

P»+, p; <&(>,.,. f , J +Q , . (4.55) 

where is the average state noise covariance for the shorter time interval, 

q = Q^±Q-. (4.56) 

Eqs. (4.55) and (4.56) are applied repeatedly until t k+1 = t c . The covariance matrix P c is a 6x6 
matrix describing the uncertainty in the Kalman filter solution at the designed time of collision. 
However, only the upper left 3x3 partition of P c is of interest, since this partition contains 


38 



uncertainties in position, and is denoted by P c . As mentioned previously, uncertainty due to 
comet outgassing is modeled as two effects. The time varying uncertainty is contained in the 
covariance Q(t ) as described above. The constant uncertainty is included by supplementing P 
with a covariance matrix calculated from differences in position at the designed time of collision 
caused by outgassing. Eq. (2.28) is numerically integrated from the time of the first observation 
to the designed time of collision, t c , once with a= 0, and then three times, each with recalculated 
using one of the o A . (i = 1 , 2, 3) previously introduced. The difference in r(t k ) between the 

a = 0 and o A . solutions are expressed as a 3x1 column matrix, 8, and the supplemented 
covariance matrix is 

p c =p t +y>,s;\ <4.57) 

1=1 

4.2.3 Position Uncertainty in the B-plane 

The covariance given by Eq. (4.57) is associated with unit vectors in the inertial coordinate 
frame, and therefore must be rotated to the B-plane coordinate frame. If C denotes the direction 
cosine matrix required for the rotation, then the covariance matrix in the B-plane coordinate 
frame is given by 

P B = CPC T . (4.58) 

The uncertainty in the orbit solution can be thought of as a 3 -dimensional error ellipsoid centered 
at the B-plane intersection point marked by the vector B. The error in position is described by a 
2-dimensional ellipse formed by the intersection of this ellipsoid with the B-plane, and is 
determined from the 2x2 partition of P B associated with the unit vectors T and R . This 2x2 
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partition is denoted by P/,. The semimajor and semiminor axes of the error ellipse, denoted by G\ 
and <x 2 respectively, are the square roots of the eigenvalues of P/,. The eigenvalues are the 
elements of a diagonal matrix D given by 


D = V r P fc V , 


(4.59) 


where V is a matrix containing the normalized eigenvectors whose directions are parallel to the 
principal axes of the error ellipse. The orientation of the error ellipse is 


9 = tan 1 


f 


V 



(4.60) 


and is measured positive clockwise from the T axis to the semimajor axis of the error ellipse. 
The out-of-plane component of the ellipsoid is associated with the error in time of arrival to the 
B -plane, denoted by 03 , and is calculated by dividing the square root of the component of Pb 
associated with S by v<*». 


4.2.4 Validation of Linearity 

Inspecting results obtained through Monte Carlo simulations will verify that the linear state 
model assumption made in the formulation of the Kalman filter algorithm is valid. Using the 
same long- and short-period comet sample orbits employed in Sec. 2.4, 100 different orbit 
solutions are obtained from 91 optical observations spaced one day apart, taken from a single 
observatory with angular resolution 77 = 0.1 arcseconds, located at the same heliocentric position 
as the Earth. The observations contain normally distributed random noise with zero mean and 
standard deviation equal to //. The samples are truncated at ±2//. Perturbations due to comet 
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outgassing are not included. Each orbit solution has a covariance matrix P&, resulting in 100 
different error ellipses centered on each of the 100 B -plane intersection points, with the 
distribution of the points themselves mimicking the shape of an ellipse. Each error ellipse is 
characterized by oi, 02 , 03 , and 6 . Average values of these parameters should agree fairly well 
with values determined from a single covariance matrix calculated from the 100 Kalman filter 
solutions of state deviation at the time of collision. Figs. 7 and 8 show the 100 B-plane 
intersection points, along with the 1 a error ellipse calculated from the single covariance matrix 
and centered on the mean coordinates of B T and B R (the 2 <7 and 3 <7 error ellipses can be 
calculated by simply multiplying the 1 (rvalues by 2 or 3 accordingly). For clarity, the error 
ellipses associated with each intersection point have been omitted, so Figs. 9 and 10 are supplied 
to show how the size and orientation of these ellipses are virtually the same for each Monte 
Carlo simulation, with only slight variations from the mean values, represented in the figures by 



Figure 7. Error ellipse for a sample long-period comet after 91 observations (la) 
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0.5 



Figure 8. Error ellipse for a sample short-period comet after 91 observations (1<7) 


a blue line. The similarity in the size and orientation of the ellipses is evidence of linearity, as 
well as the similarity of the average values to the single values of <7i, ay <T, and 6, \ listed in 
Table 1 for comparison. In other words, because the problem is linear, every orbit solution 
covariance at the designed time of collision will be virtually the same no matter where the 
B -plane intersection point is located. Since the covariance matrix can be used to describe orbit 
determination accuracy, the need to run further Monte Carlo cases is eliminated, saving time and 
computer resources. 

Since the perturbation due to outgassing is small, the problem remains linear when these 
perturbations are included. Table 2 lists values of the error ellipse parameters for the outgassing 
case. By comparing the values in Table 2 with the values in Table 1, the effect of comet 
outgassing is an increase in the size of the la error ellipse, signaling slightly less accurate orbit 
determination. 
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Figure 9. Variation in error ellipse parameters for a sample long-period comet after 91 observations 
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Figure 10. Variation in error ellipse parameters for a sample short-period comet after 91 observations 
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Table 1. Error Ellipse Parameters for Sample Long- and Short-Period Comets 



Long-| 

period 

Short- 

period 


Single 

Average 

Single 

Average 

0i (Re) 

0.914713 

0.785722 

1.158412 

0.983721 

02 (Re) 

0.224407 

0.232689 

0.188408 

0.199222 

a 3 (minutes) 

2.491799 

2.541044 

2.439856 

2.313072 

G (degrees) 

11.334516 

9.221617 

9.591973 

8.860158 

B t 


0.041637 


0.030882 

b r 


0.028503 


0.030710 


Table 2. Error Ellipse Parameters for Sample Long- and Short-Period Comets 

with Outgassing Included 



Long-| 

period 

Short- 

period 


Single 

Average 

Single 

Average 

0i (Re) 

0.923522 

0.796227 

1.166215 

0.992968 

02 (Re) 

0.229514 

0.236704 

0.193971 

0.204095 

a 3 (minutes) 

2.537502 

2.585841 

2.490257 

2.366148 

G (degrees) 

10.903420 

8.676846 

9.329242 

8.504838 

B t 


-0.068534 


-0.084936 

b r 


0.014652 


0.016525 


4.2.5 Collision Probability 

With the proof presented in the previous section that the linearity assumption is valid, Monte 
Carlo simulations are no longer required. The nominal comet trajectory is used as the reference 
trajectory in the Kalman filter orbit determination sequence, and as such, passes through the 
center of the Earth. Under this formulation the metric for describing how well the orbit has been 
determined is not an actual miss distance, but rather an Earth-collision probability. The 
projection of the Earth onto the B-plane is a disk, and the probability of impact can be computed 
by integrating the bivariate Gaussian probability density function, 
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over the effective cross-sectional area of the Earth, which accounts for the gravity of the planet 
that was neglected in the construction of the B-plane. The effective collision radius is given by 
Eq. (8.3-30) in Ref. [5] as 



where Re is the physical radius of the Earth and jue is the gravitational parameter of the Earth. 
The variables o T and o R in Eq. (4.61) are the standard deviations in B T and B R and are related to 
G\ and 02 through the expressions 

a T = tJ (cTj cos 0) 2 + [a 2 sin 9 ) 2 , 
a R = -yj (<J, sin 0) 2 + (cr 2 cos 0)~ . 

The variable p RR is the correlation coefficient of Oj and o R . 

Numerical integration of Eq. (4.61) is performed by constructing a coarse grid on the 
B-plane. At each course grid point, Eq. (4.61) is evaluated at every point on a fine grid that 
approximates a circle of radius R E , and is multiplied by the area of the fine grid square. The 
probability that the comet passes through the B-plane somewhere inside this circle is simply the 
sum of these products. Contours of constant f TR are ovals, as in Fig. 11, which shows contours of 


(4.63) 

(4.64) 
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constant collision probability for the long-period comet sample case without outgassing. The 
maximum value of collision probability occurs at the origin, and the figure can be interpreted as 
follows: if the calculated orbit crosses the B-plane at the origin (signifying a predicted collision 
at the center of the Earth), the probability of collision is 0.796; if the calculated orbit crosses the 
B- plane at (2,0) where it intersects the contour labeled -1, the probability of collision is 0.1; 
points outside the contour labeled -6 signify less than a one- in-a-million chance that a collision 
will occur. Therefore for this sample case, with 91 optical observations taken from a single 
observatory once per day, the confidence in predicting an Earth collision can never be greater 
than 0.796, or 79.6%. When perturbations due to comet outgassing are included, this probability 
decreases to 78.9%. 

Fig. 12 shows maximum collision probability, denoted by v 0 , as a function of tracking 

duration from a single observatory for both the long- and short-period sample comets, with and 
without perturbations due to outgassing included. This figure clearly shows the collision 

3 
2 
1 


-1 

-2 
-3 

-6 -4 -2 0 2 4 6 

b t (R e ) 

Figure 11. Earth collision probability for a sample long-period comet after 91 observations 
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probability results obtained from 15 to 100 days of observations, where Kq reaches approximately 
0.9. However, the results of interest are at longer tracking durations where the curves appear to 
level off and the collision probability is 0.99 or above. It is also difficult to distinguish any 
significant effect outgassing has on the results. Plotting 1 — K 0 versus tracking duration on a 
logarithmic scale, as in Fig. 13, makes the results much easier to interpret. It is clear from this 
figure that outgassing does indeed introduce uncertainty in the orbit solution. Fig. 14 shows the 
B-plane collision probability for the 141-day tracking interval. Comparing with Fig. 11, the 
ovals are now smaller in size, and the collision probability at the point (2,0) is now less than 10" 6 . 


Long-Period Comet 



Figure 12. Maximum collision probability as a function of tracking duration for 
sample long- and short-period comets 
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Figure 13. Maximum collision probability as a function of tracking duration on a logarithmic scale for 

sample long- and short-period comets 


3 - 


2 - 



-4 -3 -2 -1 0 1 2 3 4 

B 7< R E> 

Figure 14. Earth collision probability for a sample long-period comet after 141 observations 
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5 Observatory Location and Observation Schedule 


The results for the sample cases of the long- and short-period comets have been used thus far to 
aid in the presentation of the concepts and the method used to describe how well an Earth- 
impacting comet orbit can be determined. B-plane error ellipses and collision probabilities have 
been compared for these two cases, but general conclusions about the best observatory location 
or tracking schedule cannot be made based on these two cases alone. Therefore, a family of 
comets on Earth-impacting trajectories is created with varying inclinations to serve as a 
representative set of the real-life population of possible impactors, and orbit determination 
accuracies are analyzed for this family of comets. The orbits are determined from angular 
measurements over several different tracking intervals, and from single or multiple space-based 
observatories in circular orbits around the Sun at varying heliocentric distances. Radar 
measurements of range are also included in some cases. The results of this analysis will make it 
possible to compare various telescope configurations and observation schedules to be used for 
determining the orbits of Earth-impacting comets. 

The family of comets on which the remainder of the analysis will be focused contains 
five orbits with inclinations of 5°, 45°, 90°, 135°, and 175°, chosen to characterize low and high 
inclination orbits and direct and retrograde motion. Inclinations of 0° and 180° are not included 
because these cases will cause all of the elements in the first, fourth, and fifth columns of the 6x6 
matrix on the left hand side of Eq. (4.9) to be zero, and finding a preliminary orbit solution will 
not be possible. Based on the results of the sample comets in the previous section, increasing the 
orbital period (or increasing the aphelion distance from say, 40 to 300 AT!) does not cause a 
significant change in collision probability, therefore all orbits in this representative set have an 
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aphelion distance of 300 AU and a perihelion distance of 0.1 AU. The comet orbits also have Q, 


CO, and v 0 in common, with values of 30°, 143.6°, and -164.8°, respectively. 

5.1 Single Observatory at Earth 

The most convenient location for an observatory would be on the Earth, or in low Earth orbit. In 
orbit determination simulations performed for this study, it is assumed that the observatory is 
located at the center of the Earth because the distance between the observatory and the comet 
during the tracking interval is so large (approximately 5 AU). All observations are assumed 
possible; situations where the Sun or other bodies would interfere with observing the comet are 
ignored. 

Fig. 15 contains plots of collision probability, K, as a function of tracking interval for a 
single observatory coincident with the Earth. The probability was calculated at three selected 





Figure 15. Collision probability results from a single Earth observatory 
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points in the B-plane for each of the comets in the test set: at the center of the Earth, and at the 
surface of the Earth in both the Bt and Br directions. Collision probability at the center of the 
Earth is denoted by Kq, and at the surface of the Earth along the By and Br axes, collision 
probability is denoted by k t and k r respectively. One obvious characteristic of the curves in 
Fig. 15 is that each approaches a maximimum value of K at long tracking durations. This 
behavior is explained by visualizing collision probability as the area of an ellipse contained 
inside a disk. In the limit, this area approaches 1 if the ellipse center is located at the center of 
the disk, as shown in the plot on the left of Fig. 16. The area of the same ellipse with its center 
located at the surface of the disk approaches 0.5, as shown in the middle plot (half of the ellipse 
lies inside the Earth disk, and half lies outside). However, if the same ellipse is located at some 
other point on the surface of the disk as shown in the plot on the right, the area inside could be 
less than 0.5, depending on the size and orientation of the ellipse. In the limit the ellipse would 
become sufficiently small so that the area inside would approach 0.5, no matter where on the 
surface of the disk it was located. Referring back to Fig. 15, Kq approaches the limit at 1, and 
Kt approaches the limit at 0.5. The limit is approached much more slowly for Kr due to the 
orientation of the error ellipse. It is difficult to discern the variation in /rat these limits, so the 




Figure 16. Collision probability for central and limb impacts 
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behavior is shown in detail by plotting 1 -Kq, 0.5 -k t , or 0.5- on a logarithmic scale, as in 
Fig. 17. If collision probabilities of 0.99 are of interest, then approximately 111 observations 
from a single telescope at 1 AU are needed to accurately determine the direct comet orbits in the 
test set. For the retrograde orbits, 130 to 140 observations are needed to achieve the same level 
of accuracy. This difficulty is due to poor orbit geometry between the observatory and the 
comet. Fig. 18 shows the heliocentric orbits of the comet (red) and the observatory (blue) over a 
121-day tracking interval for the direct comet orbit with i = 5° and the retrograde comet with 
i = 175°. The geometry for the direct comet case results in larger heliocentric parallax between 
the observatory and the Sun-comet line than for the retrograde comet case. This larger parallax 
is beneficial to orbit determination accuracy because it helps determine the distance between 
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Figure 17. Collision probability results from a single Earth observatory on a logarithmic scale 
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( = 5° 


i = 175° 


c 

v AU 

Figure 18. Orbit geometry between an observatory at 1 AU and direct and retrograde comets 

the comet and the observatory along the line of sight. The subject of parallax will be discussed 
in more detail in the next section. 

As previously explained, Kj and Kr are sensitive to the orientation of the error ellipse. As 
the tracking duration lengthens, the orientation of the error ellipse changes so that the angle 
between the semimajor axis and the direction of the heliocentric velocity of the comet increases. 
This change in orientation is the reason why the k t curves in Fig. 17 have a “bump”. With a 
small number of observations, the majority of the orbit determination error is along track (in a 
direction parallel to the heliocentric velocity of the comet). As more observations are collected 
the error in this direction decreases, and the orientation of the error ellipse changes. Fig. 19 
shows this behavior. On the top is the error ellipse for the 45° inclination comet after 61 
observations, and on the bottom is the error ellipse for the same comet after 121 observations. 
Also shown is the direction of the heliocentric velocity of the comet projected into the B-plane, 
represented by the green lines. The error ellipse on the top is much larger and the semimajor 
axis is oriented along the direction of the velocity of the comet, signaling that the majority of the 
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Figure 19. Orientation of error ellipses for a long-period comet with 45° inclination 



error is along track. On the bottom, the error ellipse is smaller, and the angle between the 
semimajor axis and the velocity is larger, evidence that as the tracking duration increases the 
importance of the error along track decreases. 


5.2 Single Observatory in a Circular Orbit at 1 All 

Although an observatory located at the center of the Earth approximates a convenient location, 
such as low Earth orbit, there may be other positions along the heliocentric orbit of the Earth that 
result in better orbit determination accuracy. One way to investigate this is to examine the size 
of the error ellipse resulting from observatory locations that vary in increments of 5° around the 
heliocentric orbit of the Earth. Tracking durations of five days are used to determine the orbit 
solution and the associated error ellipse for each 5° increment. The results are presented in 
Fig. 20, which shows the semimajor axis of the error ellipse, <j\, versus the angle, % in the 
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ecliptic plane between the observatory and comet position vectors at the time of the first 
observation. Smaller values of 0 \ represent orbit solutions with less uncertainty than those with 
larger values of Oy. There is similarity in the variation of 0 \ between the i - 5° and i = 175° 
cases, and also between the i = 45° and i = 135° cases. Also, the comet with i = 90° has the least 
amount of variation in o\ due to the fact that the orbit plane is perpendicular to the ecliptic plane, 
and the comet appears to be directly overhead from the point of view of the observatory. 
Minimum values of G\ occur when yis approximately 0° (or 360°), with a secondary minimum 
near y= 180°; peaks in <j\ occur when yis near 90° or 270°. These variations in <j\ are due to the 
time rate of change of heliocentric parallax, d-fi/dt , given by 




/= 90 ° / = 1 35 ° / = 1 75 ° 





Figure 20. Semimajor axis of error ellipse after 5 observations as a function of the angle between the 

observatory and the Sun-comet line 


55 


Minimum values of <Ti occur at maximum values of cl/3 /dt , where yis near 0° or 180°, as shown 
in Fig. 21. Points of minimum and maximum dp/dt are labeled with the corresponding value 
of y At the points of maximum parallax rate, the observatory is moving in a direction 
perpendicular to the Sun-comet line. The maximum values of dP/dt indicate that these 
locations produce the most apparent motion between the observatory and the comet. Because the 
distance between the observatory and the comet is so large during the tracking interval, the 
comet appears to be stationary. Therefore, the main contributor to the change in parallax is the 
motion of the observatory. 



/ = 90 ° 


/= 135 ° 


/= 175 ° 




Figure 21. Semimajor axis of error ellipse after 5 observations as a function of time rate of change of parallax 


The results presented in Fig. 21 are confirmed by studying parallax rate through an analytical 
approach, where the comet is placed on the i x axis at a distance of 5 AU and dfi/dt is 
calculated at 5° increments around circular, heliocentric observatory orbits of various radii. The 
plot on the left of Fig. 22 shows the relationship between dfi/dt and y for these hypothetical 

cases. Maximum \df}/dt\ occurs when y= 0°, with a secondary maximum when y= 180°; 
minimum \dfi/dt\ occurs when y is between ±45° and ±90°. In particular, the curve 
corresponding to an observatory radius of 1 AU is at maximum \dfi/dt\ when y= 0° or 180°, and 
minimum \dfi/dt\ when y is approximately 80° or 280°. These results agree well with those 



Figure 22. Time rate of change of parallax as a function of observatory orbit radius for a hypothetical comet 

with i = 0° 
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presented in Fig. 21. Moving the observatory orbit inside 1 AU to 0.5 AU, or outside 1 AU to 3 
AU, increases the maximum parallax rate, as shown in the plot on the right of Fig. 22. Since 
better orbit determination is obtained when dfi/dt is maximum, the best locations for an 
observatory are in an orbit with a small radius, where the observatory would be traveling very 
fast around the Sun, or in a larger orbit closer to the comet. However, the cost of delivering an 
observatory to such orbits may not be justifiable based on the small increases in d/3/dt shown in 
Fig. 22. An observatory at 1 AU could perform just as well if the tracking duration was 
sufficiently long, if additional observatories were available to collect observations, or if the 
observations also included measurements of range. 

Fig. 23 shows maximum collision probability, ^ versus tracking duration for four 
different observatory orbits, corresponding to Mercury, Venus, Earth, and Mars orbit radii. The 
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Figure 23. Collision probability as a function of tracking duration for Mercury, Venus, Earth, and Mars 

observatories 
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results for short tracking intervals, when the main factor in orbit determination accuracy is 
parallax rate, agree with the conclusions drawn in the previous discussion. In each case the 
Mercury and Venus observatories out perform the Earth observatory, and the Earth observatory 
out performs the Mars observatory. However, as the tracking duration increases the orbit 
determination accuracy is improved by the additional observations, and the Earth observatory 
performs just as well, or better, than the others for direct comet orbits, as shown in Fig. 24. For 
retrograde comet orbits the comet is moving in a direction opposite to that of the observatory, 
making the parallax rate larger, and the Mercury and Venus observatories still produce the best 
orbit determination for a majority of the tracking interval for these cases. However, after 
approximately 170 days the Earth observatory results in the same orbit determination accuracy as 
the Mercury and Venus observatories. In all cases the Mars observatory produces the worst orbit 
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Figure 24. Collision probability as a function of tracking duration for Mercury, Venus, Earth, and Mars 

observatories on a logarithmic scale 
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determination accuracy, which is consistent with the hypothetical results presented previously. 
The orbit radius of Mars is near the minimum of the curve on the right of Fig. 22. Since the 
Earth observatory can be made to achieve the same orbit determination accuracy as observatories 
with higher parallax rates, there is no need to expend the additional cost of transferring an 
observatory to an orbit outside a heliocentric radius of 1 AU. 


5.3 Multiple Observatories in Circular Orbits at 1 AU 

The direction from which the comets are traveling on their Earth-impacting trajectories will not 
be known ahead of time, and it may even be necessary to track several targets at the same time. 
It is impossible to ensure that a single observatory could be placed in a location that allows 
maximum time rate of change of parallax between it and a comet. However, the odds of 
achieving the desired geometry are increased by adding a second or third observatory in the same 
heliocentric orbit. Using a multi-observatory configuration results in simultaneous angular 
measurements, which is equivalent to measurements with large parallax and no time lapse. Also, 
the number of available observations can be doubled or tripled, which will also improve orbit 
determination accuracy. Several different observatory configurations are investigated, with 
observatories placed at the Sun-Earth Lagrange points L3, L4, or L5 shown in Fig. 25 . The L3 
Lagrange point is 1 80° ahead of the Earth, the L 4 point is 60° ahead of the Earth, and the L 5 point 
is 60° behind the Earth. 

Fig. 26 shows maximum collision probability versus tracking duration on a logarithmic 
scale for various observatory configurations. The curves corresponding to the single Earth 
observatory are the same curves shown in Figs. 23 and 24. There are three curves corresponding 
to systems consisting of two observatories, each with one observatory located at the center of the 
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Figure 25. Location of Sun-Earth Lagrange points 
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Figure 26. 


Collision probability results from observatories in circular heliocentric orbits at 1 AU 
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Earth, and the second at a Sun-Earth Lagrange point. A curve is shown corresponding to a 
configuration with three observatories, one at the center of the Earth and one each at the 
Lagrange points L4 and L5. Also shown is a single observatory located at the center of the Earth 
that is capable of collecting range measurements good to 1000 km, as well as optical 
measurements. The single Earth observatory gives the poorest orbit determination in terms of 
the time required to achieve a specified value of collision probability. In the direct comet orbit 
cases, the Earth + L4 configuration performs better than the Earth + L5 configuration. The initial 
orbit geometry is such that the second observatory at L4 is moving across the Sun-comet line, 
which increases the time rate of change of parallax. An observatory at L5 is moving in the same 
general direction as the comet and does not contribute the same beneficial increase in parallax 
rate. In the retrograde comet orbit cases, the Earth + L4 and the Earth + L5 configurations exhibit 
similar performance because both systems are moving in the same general direction as the 
comet, except during longer tracking durations when the L4 observatory is moving in a direction 
perpendicular to the Sun-comet line. The system of three observatories, Earth + L4 + L5, is 
better than both the Earth + L4 and Earth + L 5 systems, but the best configurations are the two 
observatory system of Earth + L3, and the single Earth observatory capable of range 
measurements. These systems yield very similar results for retrograde comet orbits because in 
these cases the orbit geometry is ideal. The line joining the Earth + L3 observatories is 
perpendicular to the Sun-comet line, and the parallax between the observatories and the comet is 
maximum. This configuration is just as effective as a single observatory with range 
measurement capability, and orbit solutions with very high collision probability will be 
determined within three months of tracking. For direct comet orbits, the Earth + L3 observatory 
geometry is less than ideal. In fact, a tracking duration of approximately 71 days results in zero 
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parallax. Without the parallax benefit, the Earth + L 3 configuration produces less accurately 
determined orbits than the range observatory. Based on these results, the single Earth 
observatory with range capability is the best choice for observing all the comets in the test set, 
producing orbit solutions with maximum collision probability of 0.99 with 51 days of tracking. 

As mentioned earlier, with short tracking arcs of angular measurements, the major error 
at the time of impact is along the comet velocity, i.e. along track. This error results from the fact 
that the range to the comet is poorly determined by angular measurements. The advantage of 
range data is to immediately resolve this uncertainty. An approximate relation between range 
accuracy derived from angular data can be developed based on Fig. 27. Consider two 
observatories in the same heliocentric orbit that measure the angle y to locate the comet 
assuming the comet is along the normal to the line between the observatories. The heliocentric 
distance to the comet is easily derived from the distance D, which is given by 

D = dtany. (5.2) 

Assuming both observatories have an error in y of rjy, the error in the inferred distance D is 


observatory 



Figure 27. Orbit geometry for range measurements 


63 


(5.3) 


r/ D = 7 Y d sec- y 


= V r d 


1 + 


D_ 

yd y 


This relationship is plotted in Fig. 28 for 77^ = 0.1 arcsec. If the two observatories were at the 
Earth and L3 (d = 1 AU) the range to a comet at 5 AU can be determined to about 10,000 km. So 
independent range measurements with measurement accuracy much less than 10,000 km will 
significantly improve the orbit determination. From this figure it is seen that the comet would 
have to be closer than 3 AU before the angular measurements would provide the 1000 km 
accuracy assumed for the range measurements. 



Figure 28. Error in heliocentric distance of comet 
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The results presented in this section reveal an important relationship between parallax 
and orbit determination accuracy. In all previous analysis, the distance between the comet and 
the observatory is large enough so that the comet appears to be stationary during the tracking 
duration. When the tracking time scales are small compared to the time it takes the comet to 
reach the Earth, the time rate of change of parallax is an important parameter and observatory 
orbits with smaller radii perform well. On the other hand, over longer time scales the total 
parallax is the important parameter. In this case observatory orbits with larger radii are better. 
So a trade should be made concerning the amount of time between the discovery of the comet 
and the time when a collision is predicted to occur. If observing time on the order of years is 
available, an observatory outside a heliocentric radius of 1 AU will perform well because both 
the observatory motion and the comet motion help in the orbit determination accuracy. If not 
much time is available for collecting observations an observatory inside 1 AU will perform well 
because of the higher time rate of change of parallax. 
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6 Summary and Conclusions 


This research investigated the effect of tracking duration and observatory placement on orbit 
determination accuracy for a set of Earth-impacting, long-period comets with various orbit 
inclinations. Orbit solutions were obtained using a Kalman filter algorithm, in which an 
assumption is made that the comet orbit dynamics can be described by a linear model. Orbit 
determination accuracy results have been presented in terms of Earth collision probability. 
Perturbations due to comet outgassing were included in this analysis, and have been shown to be 
negligible until the comet reaches a heliocentric distance of 2.5 AU, which is approximately 200 
days of tracking from an initial distance of 5 AU. Perturbations in comet position are less than 
1 Earth radius at a heliocentric distance of 1 AU. Monte Carlo simulations proved that the 
Kalman filter linearity assumption is valid, even when comet outgassing is present. 

It was determined that the best observatory location is one in which the rate of parallax 
change is maximum based upon a brief study concerning the relationship between the time rate 
of change of parallax between the observatory and the comet and orbit determination accuracy. 
This occurs when the motion of the observatory is perpendicular to the Sun-comet line. Circular 
heliocentric observatory orbits with radii of 0.5 AU and 3 AU both have maximum parallax rates 
that are larger than that of an observatory orbit at 1 AU. The smallest maximum parallax rate 
within the inner solar system occurs in an orbit with a radius of approximately 1.7 AU. A 
comparison of collision probability results obtained from observatories in Mercury, Venus, 
Earth, and Mars orbits confirmed these findings, and in general the Mercury and Venus 
observatories performed better than the Earth observatory, while the Mars observatory performed 
worse in all instances. However, because longer tracking intervals increased the orbit 
determination accuracy of the Earth observatory to match or exceed the accuracy of the 
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observatory orbits with smaller radii, it was concluded that the additional cost of transfer and 
maintenance for an observatory outside 1 AU was not justifiable. 

A solution is considered accurate when collision probability is 0.99 or above for Earth 
observatory orbit determination. For a single Earth observatory with 0.1 arcsecond resolution, 
four months of daily observations is required to accurately determine the orbit of an outgassing 
comet. Adding an observatory positioned at the Sun-Earth Lagrange point L 4 reduces the 
necessary tracking duration to three months, and including yet another observatory at the Sun- 
Earth L 5 Lagrange point reduces the tracking duration to approximately two and a half months. 
The best orbit determination results from a single Earth observatory with range measurement 
capability good to 1000 km, which can achieve 0.99 collision probability with only 51 days of 
tracking. Similar accuracy is obtained by a system of two observatories with an angular 
separation of 180°. 

All of the observatory configurations studied are capable of producing collision 
probabilities of 0.99 or better. The reason why one configuration is deemed “better” than 
another is because the required orbit determination accuracy is obtained in less time. In fact, all 
configurations are capable of producing collision probabilities as high as 0.999999 within the 
assumptions of the study if the tracking duration is long enough. No matter how many 
observatories are used or where in the inner solar system they are located, increasing the tracking 
interval will always improve orbit determination accuracy. 
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7 Future Work 


Confirming such high collision probability as 0.999999 will require a much more detailed 
analysis than presented here, where additional force and observation error models should be 
included. The force models should include solar radiation pressure where uncertainty in comet 
mass is a consideration. The error models should include ephemeris errors, especially if the 
comet makes a close approach to Jupiter or Mars on the inbound leg of the impacting trajectory, 
and correlation between observations due to atmospheric effects should be included. Instrument 
biases should also be included, because actual error distributions will be much more complex 
than the simple normal distribution assumed in this analysis. Finally, a more consistent 
outgassing error model should be implemented. A time-varying element of uncertainty due to 
comet outgassing was included in the covariance Q it) which was updated once per day during 
the tracking intervals when observations were being collected, and approximately once every 
two weeks during the propagation from the time of the last observation to the designed time of 
collision. An appropriate time interval between updates to the covariance Q(t) should be 
determined and applied consistently throughout the analysis. 
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Appendix 


The set of six classic Keplerian orbital elements defined in Ref. [5] are 

1. a, semimajor axis; 

2. e, eccentricity; 

3. i, inclination; 

4. Q, longitude of ascending node; 

5. co, argument of periapsis; 

6. t, time of periapsis passage. 

The first two orbital elements determine the size and shape of the orbit. The next three elements 
are Euler angles that determine the orbit orientation in three-dimensional space. The sixth orbital 
element is used to determine the position of the object at a particular time. 

If the position, r, and velocity, v, of the object are known, calculating the orbital elements 
is a straightforward process. To accomplish this, determine three important vectors: the angular 
momentum vector h, the node vector n, and the eccentricity vector e. These vectors are shown in 
Fig. Al, which shows the orientation of the orbit plane in 3-dimensional space. The angular 



Figure Al. Orbital plane orientation in 3-dimensional space 
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momentum vector is normal to the orbit plane and is given by 


h = rxv. (Al) 

The node vector is defined as 

n = i z x h , ( A2) 


and lies in both the orbit and ecliptic plane. The intersection of these planes is referred to as the 
line of nodes, and n points in the direction of the ascending node. If h is parallel to i _ , then 
n = i x . The eccentricity vector is expressed as Eq. (3.14) of Ref. [4], 


e = — (vxh)- — , (A3) 

jU r 

where r is the magnitude of the position of the object. For elliptic orbits, the vector e originates 
at the focus of the orbit and lies along a line called the line of apsides. There are two important 
points in the orbit which lie on the line of apsides: periapsis and apoapsis. The eccentricity 
vector points toward periapsis, the point on the orbit where the object is closest to the primary 
body. The distance to this point is denoted r p . Apoapsis is the point where the object is farthest 
from the primary body, and the distance to this point is denoted r a . Expressions for periapsis and 
apoapsis are given by 


and 


P 

l + e 


(A4) 



(A5) 
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where p is the semi-latus rectum. 


P 



(A6) 


When the Sun is the primary body, r p and r a are referred to as perihelion and aphelion, 
respectively. 

With the three vectors h, n, and e determined, the orbital elements are simple to obtain as 
follows: 

1. Calculate the magnitude of the eccentricity vector, i.e, 

e = (e ■ e) 1/2 . (A7) 


2. Determine the semimajor axis using the expression 




where v is the magnitude of the velocity of the object. 

3. The inclination is the angle between i. and h, therefore 


cosz = 


_ h C 


h 


(A8) 


(A9) 


where h is the magnitude of h, and 0° < i < 180° . 

4. The ascending node is the angle between i x and n, hence 


sin £2 = 



n 


cos £2 = 


n ■ z. 


n 


(A10) 
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where n is the magnitude of n. 

5. The argument of periapsis is the angle between n and e, therefore 


sinry = 


e 

e 


hxn^ 

v hn / 


n ■ e 

COS (0 = . 

ne 


6. Determine the time of periapsis passage from Kepler’s equation, 

M = n(t - t) = E - e sin E , 


(All) 


(A 12) 


where n is the mean motion given by 


n = 



(A 13) 


and M is the mean anomaly, which describes an angle that evolves linearly with time. 
The mean anomaly has no significance other than convenience of notation. The 
eccentric anomaly, E , is expressed in terms of the true anomaly, v, as 


„ c + cosv . „ Vl-e 2 sinv 

cos E = , sin E , 

1 + ecosv l + ecosv 


(A 14) 


where 


e ■ r 

cosv = , 

er 



(A 15) 


and v r is the velocity along the radial direction, i.e, 


v„ = ■ 


r ■ v 


(A 16) 
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Conversely, the position and velocity of the object can be determined from a set of orbital 
elements and a specified time. First solve Kepler’s equation for E and then solve Eq. (A 15) for 
V. Calculate the magnitude of the position, r, using the equation of the orbit, 


1 + ecosv 


(A 17) 


Next, determine the position vector using the 3-1-3 Euler angle rotation [Q, i, y/\ starting with 
I >, 0, 0], i.e., 


r = r 


(cos cos El - sin sin Q cos z) i x 
(cos^sinfl + sin^'cosQcosz) i 
(sin y/ sin z) z. 


(A 18) 


where y/- OJ + v. Calculate the magnitude of the velocity from the vis-viva integral, 



\r 


l ) 

a j 


( A20 ) 


Finally, determine the velocity vector using the 3-1-3 Euler angle rotation [£2, z, y/\ starting with 

[Vr, r ~T“, 0], i-e., 
dt 


v 


(cos y/ cos Q. - sin y/ sin El cos z ) 
(cos y/ sin El + sin y/ cos El cos z ) 
(sin y/ sin z) z. 


dy r 
+ r 

dt 


(- sin y / cos El - cos y/ sin El cos z) i x 
(-sin^sinfl + cos^cosQcosz) z 
(cos y/ sin z) i_ 


(A21) 


where v r is the component of velocity along the radial direction, 


dy/ 
r — — 
dt 


— , and h = -JJip . 
r 
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